function [rmin,rmax] = numericalBounds(restr,irs,phi,opt,optimOptions)
% Function uses a constrained optimisation routine to find the upper and
% lower bounds of the identified set of the IRFs. 
% Inputs are:
% - restr: structure representing restrictions
% - irs: structure containing IRs and possibly long-run (cumulative) IRs
% - phi: structure containing reduced-form VAR parameters
% - opt: structure containing model information and options
% - optimOptions: structure containing optimisation option

signRestr = restr.signRestr;
F = restr.F;
f = restr.f;
S = restr.S;
s = restr.s;
vma = irs.vmaGK;
Sigmatr = phi.Sigmatr;
Sigmatrinv = phi.Sigmatrinv;
H = opt.H;
ivar = opt.ivar;
jshock = opt.jshock;
signNorm = opt.signNorm;

N = size(Sigmatrinv,1); % No. of variables in VAR
mZ = size(F,1); % No. of equality restrictions
mS = size(signRestr,1); % No. of sign restrictions (excl. normalisation)

% Draw opt.Qdraws initial values of Q for optimisation.
Q0 = drawQ0(restr,phi,opt);

% Set equality constraint for optimisation (Aeq*Q(:) = beq).
Aeq = zeros(mZ,N^2);
beq = zeros(mZ,1);
restrCount = 0;

for ii = 1:N % For each column of Q

    for jj = 1:f(ii) % For each restriction on ith column of Q

        restrCount = restrCount + 1;
        Aeq(restrCount,(ii-1)*N+1:ii*N) = F(restrCount,:);

    end

end

% Set inequality constraints for optimisation (Aineq*Q(:) <= b).
Aineq = zeros(mS+N,N^2);
bineq = zeros(mS+N,1);
restrCount = 0;

for ii = 1:N % For each column of Q

    for jj = 1:s(ii) % For each restriction on ith column of Q

        restrCount = restrCount + 1;
        Aineq(restrCount,(ii-1)*N+1:ii*N) = -S(restrCount,:);

    end

end

for ii = 1:N

    if signNorm == 0 % Normalisation on A0

        Aineq(mS+ii,(ii-1)*N+1:ii*N) = -Sigmatrinv(:,ii)';

    else % Normalisation on A0^(-1)

        Aineq(mS+ii,(ii-1)*N+ii) = -Sigmatr(ii,:);

    end

end

rmin = zeros(H+1,length(ivar),size(Q0,3));
rmax = rmin;
LB = [];
UB = [];

for kk = 1:size(Q0,3) % For each initial value of Q

    QInit = Q0(:,:,kk);
    QInit = QInit(:); % Vectorise initial value of Q
    
    for hh = 1:H+1 % For each horizon

        for ii = 1:length(ivar) % For each variable of interest

            % Minimise IR subject to sign and equality restrictions.
            [~,rmin(hh,ii,kk)] = fmincon(@(Q) genIRF(Q,jshock,...
                vma(ivar(ii),:,hh),N,1),QInit,Aineq,bineq,Aeq,beq,LB,UB,...
                @(Q) Qcon(Q,N),optimOptions);
            % Maximise IR subject to sign and equality restrictions.
            % Do this by minimising negative of IR.
            [~,rmax(hh,ii,kk)] = fmincon(@(Q) genIRF(Q,jshock,...
                vma(ivar(ii),:,hh),N,-1),QInit,Aineq,bineq,Aeq,beq,LB,UB,...
                @(Q) Qcon(Q,N),optimOptions);

        end

    end
    
end

rmax = -rmax;

% Take optima over different initial values.
rmin = min(rmin,[],3);
rmax = max(rmax,[],3);

end